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Abstract. Fast-slow systems are studied usually by "geometrical dissec- 
tion" [4]. The fast dynamics exhibit attractors which may bifurcate under 
the influence of the slow dynamics which is seen as a parameter of the fast 
dynamics. A generic solution comes close to a connected component of the 
stable invariant sets of the fast dynamics. As the slow dynamics evolves, this 
attractor may lose its stability and the solution eventually reaches quickly 
another connected component of attractors of the fast dynamics and the pro- 
cess may repeat. This scenario explains quite well relaxation oscillations and 
more complicated oscillations like bursting. More recently, in relation both 
with theory of dynamical systems [11] and with applications to physiology 
[10, 26], a new interest has emerged in canard cycles. These orbits share the 
property that they remain for a while close to an unstable invariant set (cither 
singular set or periodic orbits of the fast dynamics). Although canards were 
first discovered when the transition points are folds, in this article, we focus 
on the case where one or several transition points or "jumps" are instead 
transcritical. We present several new surprising effects like the "amplifica- 
tion of canards" or the "exceptionally fast recovery" on both (l+l)-systems 
and (2+l)-systems associated with tritrophic food chain dynamics. Finally, 
we also mention their possible relevance to the notion of resilience which has 
been coined out in ecology [19, 22, 23]. 

Introduction 

Systems are often complex because their evolution involves different time 
scales. Purpose of this article is to present several phenomena which can be 
observed numerically and analyzed mathematically via bifurcation theory. 

A first approximation for the time evolution of fast-slow dynamics is often seen 
as follows. A generic orbit quickly reaches the vicinity of an attractive invariant 
set of the fast dynamics. It evolves then slowly close to this attractive part 
until, under the influence of the slow dynamics, this attractive part bifurcates 
into a repulsive one. Then, the generic orbit quickly reaches the vicinity of 
another attractive invariant set until it also loses its stability. This approach 
is a quite meaningful approximation because it explains many phenomena like 
hysteresis cycles, relaxation oscillations, bursting oscillations [13, 27, 28] and 
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more complicated alternation of pulsatile and surge patterns of coupled GnRH 
neurons [6, 7]. 

Discovered by E. Benoit, J.-L. Callot, M. and F. Diener (see [3]), canards 
were first observed in the van der Pol system: 

(o.i) e± = y - f( x ) = y - (t + %2 ) 

y = x — c{e) 
where c(e) ranges between some bounds: 

(0.2) c + exp (-^j < c(e) < c + 

A canard is an orbit which remains for a while in a small neighborhood of 
a repulsive branch of the critical manifold y = f(x) i.e. a connected set of 
repulsive points for the fast dynamics. In the following, we consider the canard 
phenomenon in its broader sense of delay to the bifurcation of the underlying 
fast dynamics under the influence of the slow dynamics. More recently, F. 
Dumortier and R. Roussarie [11] contributed to the analysis of such orbits by 
blowing-up techniques and introduced the notion of canard cycle. There are now 
several evidences showing the relevance of this notion to explain experimental 
facts observed in physiology (see [10, 26]). 

The systems presented here display anomalous long delay to ejection from 
the repulsive part of the fast dynamics. 

1. Enhanced delay and canard cycles of planar (1+1)-dynamics 
1.1. Dynamical Transcritical Bifurcation. 

The classical transcritical bifurcation occurs when the parameter A in the 
equation: 

(1.1) x = —Ax + x 2 

crosses A = 0. Equation (1.1) displays two equilibria, x = and x — A. For 
A > 0, x = is stable and x = A is unstable. After the bifurcation, A < 0, x = 
is stable and x = A is unstable. The two axis have "exchanged" their stability. 

The terminology "Dynamical Bifurcation" (due to R. Thorn) refers to the sit- 
uation where the bifurcation parameter is replaced by a slowly varying variable. 
In the case of the transcritical bifurcation, this yields: 

where e is assumed to be small. This yields: 

x = -{-et + y )x + x 2 , (yo = 2/(0)) 
which is an integrable equation of Bernoulli type. Its solution is: 

= W-rWl , {xo = x{0)) 

1 — xo Jo exp[—Y(u)\au 
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where: 

ft ft t 2 

Y (t) = / y(s)ds = / (-E8 + y )ds = -e— + y t 
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Figure 1. Orbits of system (1.2) starting from (0.05z, 1), for 
i — 1, .., 6, with e — 0.1. The critical set - formed by the singular 
points of the fast dynamics - is shown in green and red. On this 
set, the red points are repulsive for the fast dynamics and the 
green points are attractive. Each orbit reaches a neighborhood 
of the attractive manifold x = 0, y > 0, goes down along x — 0. 
When y becomes negative, despite the repulsiveness of x = 0, y < 
0, each orbit remains for a very long time close to x = before 
drifting away. 

If we fix an initial data (x ,y ), y > 0, < x < yo/2, and we consider the 
solution starting from this initial data, we find easily that it takes time t = yo/s 
to reach the axis y = 0. If xq is quite small, that means the orbit stays closer 
and closer of the attractive part of the critical manifold until it reaches the axis 
x = y and then coordinate x starts increasing. But now consider time cy^/e, 
1 < c < 2. Then, a straightforward computation shows that: 

and: 
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This shows that, despite the repulsiveness of the half-line x = 0, y < for the 
fast dynamics, the orbit of (1.2) remains for a very long time close to x = 0, 
indeed x(cyo/e) < 2xq for e small enough (see Figure 1). Note that, afterwards 
for larger time, the orbit blows away from this repulsive axis. This phenomenon, 
although quite simply explained, is of the same nature as the delay to bifurcation 
discovered for the dynamical Hopf bifurcation, see for instance [1, 5, 12, 29]. This 
well-known effect is instrumental in the systems we study in this article. Some 
related work has been done in computing entry-exit relation for the passage 
near single turning points (see [2, 9]). 

1.2. Enhanced delay to bifurcation. 

Consider the system: 

( 13 N X = (1-X 2 )(X-1/) 

^ ' ' y = ex 

The critical set - defined as the set of singular points of the fast dynamics - 
is the union of the straight lines x — — 1, x — 1, and y = x. In the following, 
we call "slow manifolds" the two straight lines x — — 1 and x — 1, as they are 
invariant for the critical system {(1 — x 2 )(x — y) = 0, y = x}. A quick analysis 
shows that, as the slow variable y (considered as a parameter) evolves, the fast 
system undergoes two transcritical bifurcations: at x — — 1 for y = — 1, at x = 1 
for y — 1. 

As we recalled in Subsection 1.1, a typical orbit starting from an initial data 
( x o,yo)i \ x o\ < 1; close to (x = —l,y > —1), first goes down along x = — 1 and 
displays a "delay" along the repulsive part of the slow manifold (x = —l,y < 
— 1). Then, under the influence of the fast dynamics, it quickly reaches the 
attractive part (x = l,y < 1) and moves upward to the other transcritical 
bifurcation point. There, it again displays a delay along the repulsive part (x = 
l,y > 1). Then, it quickly reaches (x — l,y > —1) and starts again. Hence, 
there is a mechanism of successive enhancements of the delay to bifurcation 
after each oscillation generated by the hysteresis. We proved in [1 1] the: 

Theorem 1. For all initial data inside the strip —1 < x < 1, for all 5 and for 
all T , the corresponding orbit spends a time larger than T within a distance less 
than 5 to the repulsive part of the slow manifolds. 

We also proved in [11] the: 

Theorem 2. Given any initial data (xo,yo) outside the strip \x\ < 1, the cor- 
responding orbit is asymptotic to y = x. 

1.3. Structural stability of the enhanced delay. 

The theory of the structural stability of fast-slow systems remains to be found. 
In this subsection, we actually adopt a very pragmatic approach and restrict 




Figure 2. Typical orbits of (1.3) with e = 0.1 starting from a 
point near the origin in panel 1) and from various initial datas 
outside the strip \x\ < 1 in panel 2). Double arrows are added on 
the fast parts of the orbits and single arrows on the slow parts. 
Panel 1): The delay to the transcritical bifurcation undergone by 
the orbit is enhanced at each passage. Consequently, the time 
taken to escape a given neighborhood of x = — l,y < — Ion one 
hand and a given neighborhood of x = 1, y > 1 on the other hand 
is longer for each oscillation than for the preceding one, although 
these half-lines are repulsive for the fast dynamics. 
Panel 2): Each orbit is asymptotic to y = x. 
These numerical simulations were performed with order 4 Runge- 
Kutta integration scheme (absolute error value: 10~ 9 ; relative 
error value: 10~ 12 ; mean integration step: 10" 5 ). 

ourselves to numerical simulations. We perturb the system by adding an ar- 
bitrary small perturbation generated by a chaotic system, more precisely the 
Rossler system [25]. 

/ 14 n ei = (I - x 2 )(x - y) + an(t) 

^ ' ' y — x 

where n(t) is the first variable of a Rossler system: 

n = —u — v 

(1.5) u = n + au 

v = b + {n — c)v 

with a = 0.1, b = 0.1, and c = 14 with some fixed initial data. 

We choose a quite small. The perturbation does not change the behavior of 
the typical orbit for a while. It starts to oscillate between x = 1 and x — — 1: it 
remains for a long time alternatively near each of these lines and undergoes fast 
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motions after ejection. But, as the orbit approaches closely these axes (within 
the distance of a), the small chaotic perturbation becomes operating and moves 
slightly the orbit outside of the strip. After crossing one of the straight lines 
i = -1 or i = 1, the orbit moves quickly to the axis y = x as previously shown. 

4 
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Figure 3. Example of an orbit of (1.4), with e = 0.5 starting 
from (x, y) = (10~ 3 , 0) for a = 10~ 4 . Double arrows are added on 
the fast parts of the orbit and single arrows on the slow parts. 

Figure 3 displays an example of such an orbit starting from (0.001,0) for 
a = 10~ 4 . Despite the very small magnitude of the perturbation - less than 
3.10" 3 - the orbit exits the strip \x\ < 1 before the third oscillation around the 
origin. This simulation shows the loss of the oscillatory behavior as, under the 
influence of a perturbation of very small magnitude, the orbit crosses one of 
the slow manifolds as soon as x is sufficiently close to —1 or 1. We have thus 
obtained numerical evidence that the system is not structurally stable. 

It is interesting to point out the fact that the numerical simulation of (1.3) 
highlights also this structural unstability. In fact, the integration scheme itself 
- whatever the integration step and the tolerances - provides errors. They give 
rise in the long run to intrinsic approximation: after many oscillations, the orbit 
starting from \x\ < 1 passes so close to the slow manifolds that the numerical 
integration leads, sooner or later, to the approximation x = 1 or x = — lor can 
eventually cross one of these lines. 

1.4. Canard cycle for a double transcritical system. 

In this subsection, we introduce a new system inspired by the preceding (1.3): 

/ 16 ) * = U ~x 2 )(x-y) 

y = ex(y + b)(a - y) 

where a, b > 1 are parameters. The critical set is the same as the critical set 
of (1.3), formed by the three straight lines of equation x = — 1, x — 1 and 
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y = x. The dynamics displays again two transcritical bifurcations: at x — — 1 
for y — — 1, at x — 1 for y — 1. But, the novelty is in the factors (a — y) and 
(y + b) 'my which yield bounded orbits. We restrict the phase space to the 
compact set: 

K = {(x,y)\ - l<x<l,-b<y<a} 

As for (1.3), the origin is a repulsive focus of (1.6). Except the origin, the 
singular points of (1.6) in K are the summits of the rectangular boundary of 
K, (—1, —b), (1, —b), (1, a), (—1, a), and are all of saddle type. As the straight 
lines x = —1, x = 1, y = a, y = —b are invariant under the flow of (1.6), we 
deduce the stable and unstable manifolds of each saddle: 



saddle 


stable manifold 


unstable manifold 


(-1,-6) 


x = —1, y < a 


y = —b, x < 1 




y = —b, x > — 1 


x = 1, y < a 


(l,o) 


x = 1, y > —b 


y — a, x > — 1 


(-l,a) 


y = a, x < 1 


x = 1, y > —b 



Consequently, the interior of K is invariant under the flow and contains only 
a repulsive focus. The boundary of K is a graphic (cu-limit set formed by the 
union of the saddles with their separatrices) . 

Following the study on the dynamical transcritical bifurcation, one expects 
that the value of y along an orbit oscillates between —b and min(a, b + 2) if a > b 
or between — min(6, a + 2) and a if a < b. However, the exponential attraction 
of x = —1, — 1 < y < a (resp. x = 1, — b < y < 1), formed by attractive points 
of the fast dynamics, produces a delay to bifurcation and keeps the orbit near 
x = —1, —b<y<—l (resp. x — 1, 1 < y < a), formed by repulsive points of 
the fast dynamics. The delay could be great enough so that the orbit reaches 
a small neighborhood of y = 0. Hence, after the fast motion, the orbit tracks 
the other manifold x — 1, — b < y < 1 (resp. x = —1, — 1 < y < a) even closer 
than during the preceding passage. The delay obtained afterwards is then again 
enhanced, as the orbit approaches even closer the slow manifolds. Hence, an 
orbit starting near the origin, after several oscillations, reaches alternatively a 
small neighborhood of both y = a and y = —b, whatever the values of a and b. 

Figure 4 illustrates this global behavior with fine step numerical simulations 
for various sets of parameter values. As expected, the smaller e, the faster the 
orbit reaches a given neighborhood of x = — 1 (resp. x = 1). Starting near the 
origin, if e is small (for instance, equal to 0.5 for the simulations presented in 
Figure 4), a few oscillations suffice to obtain a part close to the graphic. 

Thus, the orbits of system (1.6) generate a new type of canards. On the 
contrary of the canards discovered in the solutions of the van der Pol system, all 
the orbits of system (1.6) starting from il~\{(0, 0)} display delays to bifurcation. 
Moreover, we have a very simple way to modulate these delays by choosing the 
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Figure 4. Orbits (red) of (1.6) for e = 0.5, b = 2 and respec- 
tively a = 5 in panel 1) and a = 30 in panel 2). In both cases, 
the initial data is (10~ 3 , 0). Comparison of panels 1) and 2) shows 
that, whatever the values of a and b, the orbit enters alternatively 
small neighborhood of x = — 1, x — 1, y = — b and y = a (blue 
straight lines). These numerical simulations were performed with 
order 4 Runge-Kutta integration scheme (absolute error value: 
10~ 9 ; relative error value: 10~ 12 ; mean integration step: 10~ 6 ). 



parameter values a and b. It is worth noticing that, as in system (1.3), small 
perturbations may provoke dramatic changes in the orbits. 

Now, we explore the ejection mechanism - leading to this type of oscillations - 
using the conjugacy of (1.6) with an appropriate semi-local form and the analysis 
of the transition near the saddle (x, y) = (—1, —b) (the case of (x, y) = (1, a) is 
similar). First, we translate the origin to the saddle via the change of variables 
{X — x + 1,Y — y + b} to obtain the new form of the system: 



X = X(X-Y + b-l)(2-X) 
y ' ' Y = eY{X -l)(a + b-Y) 

Note that the left point of transcritical bifurcation of the fast dynamics (x, y) = 
(-1, -1) reads now (X, Y) = (0, b - 1)._ 

Consider an initial data (X, Y) with X g]0, 1], 6—1 < Y < a+b. As explained 
above, for e small, the corresponding orbit of (1.7) reaches quickly the vicinity 
of X — 0, y > b — 1, X decreases very quickly while Y remains almost constant. 
Hence, along this fast branch of the dynamics, variable Y is given by Y + 0(e) 
until X = 0(e). Afterwards, approximation of the dynamics near X = by the 
slow motion and direct integration yield an approximation for the time needed 
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to reach Y = b — 1 : 



■Y^b-l 



In 



(a + l)Y 



(b-l)(a + b-Y) 



e^o £ (a + b) 

The leading term of the x-component Cexp(— k/e) displays: 
(1.8) C = (b- l) 2 ^(a + l) 2 ^ 



0(1) 



;i.9) 



k 



(b - 1) \n(Y) + (a + 1) ln(a + b — Y) 



> 



We now study the delay to bifurcation and consider the entry of the orbits 
coming from above Y = b — 1. Hence, we consider initial values of X which are 
exponentially small with respect to e: C exp(—k/e). 

To specify the transition induced by the flow near the saddle, we consider the 
two sections: 

(1.10) S m = {(X,6)\0<X < v } 

(1.11) ^out = {(r],Y)\0<Y <5} 

with r\ > a small fixed parameter and < 8 < b — 1. We note: 

17 = {(X,Y)\0 < X < 77,0 <F <5} 

the rectangle delimited by S out and the stable and unstable manifolds of 
the saddle (see Figure 5). 




Figure 5. The transition function induced by the flow of (1.7) 
is well-defined from Ej n into T, out . 

Hence, as Y < in {(X,Y)\0 < X < 1,0 < Y < X + b - 1} and 5 < b - 1, 
any orbit starting from E in enters U and escapes from U through T, out . Thus, 
the transition function induced by the flow is well-defined from Ej n into T, out . 
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As X is small in U (smaller than 77), system (1.7) is conjugated to: 

X = 2(b-l)X-2XY 
[ } Y = sY 

Note that (X, Y) = (0, 6 — 1) is also a point of transcritical bifurcation for this 
system. Direct integration provides the orbit from (Xq,5): 



X(t) = X exp 
Y(t) = 5e~ £t 



2(6 — 1)* 



25 



e~ £t - 1 



As explained previously, we consider initial data of type X = Cexp(— k/e), 
where k, C > 0. If r\ < C, for e small enough, i.e.: 



k 



e < 



(X , 5) lies in S in . If 77 > C, as s, k > 0, all values of e fulfill this property. 

The time T needed to go from (X ,S) to (77, Y(T)) G S oltt along the flow is 
the solution of: 

2(6 - I T + — e" eT = - + In -L + — 
£ e C e 

The solution T can be expressed via the "Lambert function" Wc - inverse 
function of w —> we w (see [8]). This yields: 



1.13) T = -W c 



5 



b-1 



exp 



e\n{ V /C) +k + 25' 



2(6 -r 



+ 



In 



77 k + 25' 

c + — 



As expected, the transition time is 0(l/e) and the leading term of T is: 



1.14) 



1 

e~*0 e 



5 



2S+k 



b-1 



e + 



25 + k 



0(1) 



This displays the transition function (Xq,5) — > (r],Y out (Xo)) with: 
1.15) Y mt {Ce- k s) = 



5 exp 



-Wc 



5 



b-1 



exp 



eln(77/C) + k + 25' 



2(6-i; 



6-1 



' v k + 25' 
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This shows that, even if the transition time tends to +00, the exit function 
displays as 0(l)-leading term: 



;i.!6) Y out (Ce 



e-*0 



£ exp 



-Wr - 



1 



_ 2<S+fc 

-e 2(b-i) 



25 + k 



0(e) 



Hence, the longer the orbit has remained near the slow manifold, the stronger 
the ejection (to reach the neighborhood of (f/, y o „ f ) given by (1.16)) is. This 
property of the orbits - staying a very long time near the repulsive part of 
the slow manifold without squashing on the y-nullcline - together with the 
increasing strength of ejection is what we call the "exceptionally fast recovery" . 
Similar study can be done for the other transition of (1.6) and this leads to 
similar results (where a and b are exchanged). 

Finally, note that, for given values of parameter a, b, e, the values of param- 
eters k and C in (1.14) and (1.16) are approximated for the global dynamics 
using (1.8) and (1.9). 



2. Saddle-node transcritical ejection in a 
prey-predator-superpredator model 

2.1. Tritrophic food chain dynamics. 

In the late seventies, interest in the mathematics of tritrophic food chain mod- 
els (composed of prey, predator and superpredator) appear (see, for instance, 
[15, 16]). Related predator-prey models with parasitic infections were studied 
later [17]. In the nineties, in [18, 24] and [20], the existence of chaotic attractors 
was discussed. There are many more recent contributions, that we can not refer 
in more details (see, for instance, [21, 28]). 

We investigate here the following: 

du u(r(i- u ^ W 



dT V V KJ Br + U 

(2A) % = yU^L t a 2 w 



dT V B 1 + U B 2 + V 

dW ( A 2 V 

eW E 2 — — - D< 



dT V B 2 + V 

which represents the interactions between three populations U, V and W. The 
variable U stands for the prey, V for its predator and W for a superpredator 
of V. The threshold constant K > and the intrinsic growth rate of the prey 
R > characterize the logistic evolution of U. 
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The predator-prey interactions are described by two Holling type II factors 
defined by the positive parameters: 

Aj : the maximum predation rates 
Bj : the half-saturation constants 
Dj : the death rates 
Ej : the efficiencies of predation 

j = 1 relates to the predator V, j = 2 relates to the superpredator W. It 
is assumed that the evolution of the superpredator is slower than those of the 
predator and the prey. Then we introduce different time scales by means of the 
constant < e < 1. 

2.2. Bifurcations of the fast dynamics. 

The global behavior of this system (existence of global periodic orbit, bifurca- 
tions of limit cycles, early crisis in the predator membership) has been studied 
in [27] and [28]. We recall briefly the classical situation of interest for us for 
which the system is bistable. 

In order to obtain a simpler and more useful analytic form, Klebanoff and 
Hastings proposed in [20] the following rescalings: 

U V w 

which yields: 

x = x (l -x - a \ ) = f(x,y,z) 

(2-3) y = y\YTT tli - 1 = 9(x,y,z) 

\l + bix 1 + b 2 y J 

* = £z ( T~^T d 2} = h(x, y, z) 

\l + b 2 y J 

where a,, bj, dj,j = 1,2 are positive parameters (see [20] or [27] for their expres- 
sions in function of Aj, Bj, Dj, Ej, R and K. All axes and faces of the positive 
octant M? + are invariant sets of (2.3). Thus, we limit the phase space to this 
positive octant. 

Considering the slow variable parameter, we describe the sequence of 

bifurcations undergone by the two-dimensional fast dynamics. To this purpose, 
it is convenient to introduce the critical set: 

(2.4) C = {{x, y, z) G R%\f(x, y, z) = 0, g{x, y, z) = 0} 

formed by the singular points of the so-called Boundary-Layer System (BLS), 
obtained from (2.3) with e — 0. For any point (x,y,z) G C, (Sc,y) is a singular 
point of the fast dynamics for z — z. In the following, we note tt the projection 
from M.+ into M^_: ir(x,y,z) = (x,y). 
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First, we assume: 
(2.5) G = ai-di(l + 6i)>0 

to ensure that the singular point: 

f d t G \ 



(2.6) 



G + di (G + dx) 



2 ' 



lays in the phase space 

The critical set writes C = A U £ where: 

(2.7) A = {(1,0, z)\zeR + } 

(2.8) £ = {(x,yc(x),z c (x)) ER 3 + \x E[0,1}} 
and: 

(2.9) y c (x) = -(l-x)(l + b lX ) 



(2.10) z c (x) 



[a\x — d\ (1 + bix)) (ai + 6 2 (1 — x)(l + Oi^)) 
0102(1 + Oix) 

Note that A always intersects £ at the point T = (1,0, Zr) where: 

(2.11) z T = jz 7-t > 

In the following, we assume that di is small enough so that: 

(2.12) there is a point xp > so that ^(xp) = 0. 

We note P = (xp,yp, zp) = (xp,yc(xp), zc(xp)). Let us remark that, under 
assumption (2.12), zp < z P . Consequently, £ is fl-shaped in (cf. Figure 6) 
and we note: 

• for each z G ]z T , zp[, S z the unique point (x, yc(x), z) G £ 
such that x > xp ; £s = U \S Z } ; 

• for each 2 G [0, Zp[, R z the unique point (x, yc(x), z) E £ 
such that x < Xp ; £± = U |-R 2 i- 

0<^<zp 

Thus, £ = TU£ s UPU/; ± (see Figure 6). 

The points 7r(S z ) are saddle type singular points of the fast dynamics, (1, 0) is 
a saddle for z < zp and an attractive node for z > Zp. Additionally, we assume 
that d\ is small enough such that there exists a point zh g]zt, zp[ so that: 

for < z < Zh, k(R z ) is a repulsive focus, 
for Zh < z < zp, ir(R z ) is an attractive focus. 

Figure 6 displays C and its splitting according to the nature of the singular 
points for the fast dynamics. It can be seen as a bifurcation diagram of the 
fast dynamics as the bifurcation parameter z varies. Hence, as z decreases, the 
following sequence of bifurcations occurs: 

• for z > zp, the attractive node (1,0) is the unique singular point. 
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Figure 6. Critical set of (2.3) (in red) - set of singular points 
of the Boundary-Layer System (BLS) obtained by setting e — 0. 
It is splitted according to the nature of the singular points for 
the fast dynamics with the corresponding value of z: A_ and 
formed by attractive nodes, and C$ by saddles, C + by 
repulsive foci. The points P, H, T are respectively the points of 
saddle-node, Hopf and saddle-node transcritical bifurcation of the 
fast dynamics. For z e [0, zh[, the repulsive focus tt(R z ) e £+ is 
surrounded by an attractive limit cycle of the fast dynamics. The 
union M. of these limit cycles (in green) is an invariant attractive 
manifold for the (BLS). 

as z = zp, an inverse saddle-node bifurcation occurs at (xp,yp). 

for zh < z < zp, there are three singular points: the attractive node 

(1,0), the saddle n(S z ), the attractive focus ir(R z ). 

as z = zh, tt(R z ) undergoes a supercritical Hopf bifurcation. 

for zp < z < zh, there are three singular points: the attractive node 

(1,0), the saddle tt(S z ), the repulsive focus n(R z ) surrounded by an 

attractive limit cycle, born from the Hopf bifurcation. 

as z — zt, a saddle-node transcritical bifurcation occurs at (1,0) {tt(S z ) 

and (0, 1) exchange their stability). 

for < z < zp, there are two singular points in the positive octant: the 
saddle (0, 1) and the repulsive focus surrounded by an attractive limit 
cycle. 
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We note M. the union in of the planar attractive limit cycles surrounding 
the points R z for all zh < z < zp. A4 is thus an attractor for the (BLS). 

A precise study of this sequence of bifurcations is available in [27, 28], includ- 
ing the double homoclinic bifurcation that may occur in a certain range of the 
parameter space. However, here, we focus on the oscillatory behavior of the sys- 
tem in the phase space and control the two types of ejections giving rise to the 
hysteresis loop: the saddle-node bifurcation and the saddle-node transcritical 
bifurcation. 

2.3. Global behavior of the system and canard cycles. 

The separatrix z = is the plane of equation: 

d 2 

a 2 - o 2 d 2 

For < d 2 < a 2 /b 2 small enough, this plane separates A near which z < and 
Ai near which i > 0. 




Figure 7. Outline of a typical orbit of system (2.3) generating 
bursting oscillations. Double arrows are added on the fast parts 
of the orbit and single arrows on the slow parts. 

Hence, any orbit of system (2.3) starting near A_ goes down along A_. Once 
z < zti it undergoes a delay to bifurcation and keeps tracking As for a while, 
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although this branch is formed by saddles of the fast dynamics. After this 
delay, the orbit quickly reaches the vicinity of Ai. It goes up while spiraling 
around A4, then around As z becomes larger than zp, the orbit reaches 
the vicinity of A_ and repeats the same sequence of motions. Figure 7 displays 
the geometric invariants of the fast dynamics and a schematic orbit of (2.3). In 
this setting, the delay to the transcritical bifurcation, that the orbit undergoes, 
qualifies the terminology of canard cycle. 

Figure 8 displays a numerical simulation of a typical orbit of (2.3) with: 

ai = 0.8, 6i = 4, di = 0.1, e = 0.1, 
a 2 = 42, b 2 = 40, d 2 = 1. 




2.5 



Figure 8. Orbit of system (2.3) with the parameter values given 
in the table of Subsection 2.3. Double arrows are added on the 
fast parts of the orbit and single arrows on the slow parts. 

2.4. Local form near the saddle. 

The most surprising behavior of these orbits is the long time needed to escape 
a neighbordhood of the point (1, 0, 0) and the "exceptionally fast recovery" of the 
variable y. In this subsection, we perform successively two changes of variables 
to obtain an appropriate local form near the saddle. 

The Jacobian matrix associated with the system (2.1) at the singular point 
(x, y,z) = (1,0,0) reads: 

H -TTfe \ 
V -ed 2 J 
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It admits — 1, 



di 



G 



and — ed 2 as eigenvalues. Under the assumption 



1+61 "' i i+fci 

(2.5), the singular point (1, 0, 0) is then a saddle with a two-dimensional stable 
manifold and a one-dimensional unstable manifold. 
Eigenvectors associated with the eigenvalues are: 



eigenvalue 


eigenvector 


-1 


(1,0,0) 


G 
I+61 


(_l,i(G + l + 6i),0) 


-ed 2 


(0,0,1) 



Set: 



a 



K 



1 



h){d x 



The parameter a is positive under the assumption (2.5) (G > 0). 
Hence, via the change of variables: 



x 

y 

z 



1-X 

aY 

Z 



Y 



the singular point (x, y,z) = (1,0,0) translates to (X,Y,Z) = (0,0,0). Now, 
the eigendirections of the saddle coincide with the (X, Y, Z)-axis. In this new 
system of coordinates, (2.3) reads: 



X 



(2.13) 



Y 



"(1 
+Y 

Y 



-X 
-1 



Y)(X + Y) 

(x-y)(d! 



a 2 Z 



1 + 61(1 ^ 

ai {l + X-Y) 



Ze 



f 6i(l - 
a 2 aY 



1 + 6 9 «y 



X-F) 
-d 2 



X-Y) 
d\ — 



1 + b 2 aY 
a 2 Z 



boaY 



This allows to consider X > in the following. 

Actually, an eigenvector associated with each saddle (0, 1) of the fast dynamics 
writes: 



(2.14) 



a 2 Z(l + bi) 



•1 



.l + fti + G-aa^l + fti, 
This vector is well-defined in the positive octant as long as 



Z < 



l + h + G 
a2(l + &0 



1 

a 2 



G 



61 



Thus, under the assumption (2.5), it is well-defined at least for Z G [0,zt] 
and (3{Z) > 0. Note that, for Z = zt, this vector gives the central direction 
associated with the eigenvalues of the non hyperbolic point (0, 0) of the fast 
dynamics. 



18 



ALEXANDRE VIDAL AND JEAN-PIERRE FRANQOISE 



For each Z e [0, Zr], the vector (2.14) gives the tangent line to the unstable 
manifold of the saddle (0, 0) of the fast dynamics with Z = Z. Thus, we have 
an approximation for X and Y small of the invariant manifold of (2.13) with 
e = 0. As this manifold is normally attractive, it persists for e > small enough 
into a normally attractive manifold of (2.13). We proceed with another change 
of variables: 

u = X-f3(Z)Y 
y = Y 
z = Z 

locally near X = Y = 0,0 < Z < zt, system (2.13) reads: 

u = —u + 0(ev,uv,u 2 ) 

(2.15) v = G'v — a 2 vw + vO(u,v) 

w = —ed 2 w + wO{v 2 ) 

where: 

G' = ^— = ^--d 1 
l + h l + h 

The perturbed attractive manifold in this new set of parameters is approximated 
by u = 0, < Z < zt. 

The flow on the invariant perturbed manifold is conjugated - locally near 
u — v = - to the system: 

(2.16) v = v( ^— di) — a 2 vw 



+ 

w = —ed 2 w 

After setting w = 2Y/a 2 and replacing ed 2 by e, we recognize the normal form 
(1.12) introduced in Subsection 1.4 with parameter values G' — (b — 1) > 0. 

2.5. Transition function: recovery of the predator. 

Consider the two sections: 

(2.17) E in = {(u,v,6)\0 <v < V } 

(2.18) E out = {(u,7 1 ,Y)\0<u<£,0<w<5} 

where 77, £ > are small but fixed parameters and < 5 < z T (see (2.11)) is 
fixed. We note: 

U = {(u, v,w)\0 < u < £, < v < 77, < w < 5} 

The preceding analysis shows that, for e small enough, an orbit of (2.16) starting 
from E in enters U and exits from U. Hence, the transition function induced by 
the flow of (2.16) is well-defined from S« n into T, out . 
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As the typical orbits of system (2.3) enter an exponentially small neighbor- 
hood of A 5 (see Figures 7, 8 and Subsection 2.3), we restrict initial data of 
(2.16) to be exponentially close to X = Y = in £ in . Thus, consider: 



(2.19) 
(2.20) 



Uq 
Vq 



C 2 e- 



By direct application of the local analysis made in Subsection 1.4, we obtain 
the leading term of the transition time from (u ,v ,S) G S in to S out along the 
flow of (2.16): 



Wr 



doS <^2 6 + k 2 d 2 

e g 



G' 



+ 



a 2 5 + k 2 d 2 
G' 



+ 0(11 



The transition function (uq,vq,S) 
fled as well: 



{u out (uo,v ),r),w ou t(u ,v )) can be speci- 



(2.21) u out {C 1 e--,C 2 e 



-Wr 



CL 2 5 a-2 s + k 2 d 2 

e G 1 

G' 



a 2 5 + k 2 d 2 



(2.22) w cmt (C 1 e-^,C 2 e-^) 
= 5ex\ 

It is worth noticing that the it-component of the transition function tends 
actually to as e — > 0, but that the w-component experiments a basal threshold, 
even if this threshold is very low. 

2.6. Interpretation in terms of ecological parameters: loss of resilience. 



-Wr 



a 2 5 o, 2 s+k 2 d 2 \ a 2 5 + k 2 d 2 



G' 



doG' 



+ 0(e) 



In contrast with system (1.6), the two ejection mechanisms in (2.3) are pro- 
vided by different types of bifurcation of the fast dynamics: a saddle-node tran- 
scritical bifurcation at T and a saddle-node bifurcation at P. Near P, canards 
could eventually appear for which the orbit tracks the manifold Cs formed by 
saddle of the fast dynamics. However, this type of canards can occur only if P is 
near the plane z = (within a 0{\/e) distance). But the assumption that this 
plane separates properly A and the manifold M. (d 2 small enough) generally 
forbids such behavior. 

Hence, any orbit in the strictly positive octant displays, sooner or later, the 
motion described in Subsection 2.3 as there exists a globally attractive periodic 
orbit of this type (see [27, 28]). Thus, the entry of this periodic orbit - as the 
entry of any orbit after spiraling around M. - in a given neighborhood of A 
occurs for z = zp + 0(e 2 ^ 3 ). The delay to bifurcation after the passage in the 
vicinity of T is then mostly characterized by zp — zt- In fact, as calculated 
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for (1.6) in (1.9), the contraction exponent between A and the orbit is chiefly 
determined by z P — z T . 

The "exceptionally fast recovery" that we have pointed out therein may have 
some connection with the notion of "loss of resilience" which has been discussed 
in several dynamical models related to Ecology [19, 22, 23]. The long return 
times associated with a loss of resilience are caused by slow dynamics near the 
unstable equilibrium. We expect that new developments of this subject could 
possibly benefit of the notion of canards and canard cycles. 

Acknowledgements: The numerical simulations were performed using XPP- 
AUT: http://www.math.pitt . edu/~bard/xpp/xpp .html 
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